module upwind_mod
	
  use constants_mod
  use namelist_mod
  use boundary_mod, only: pbc_idx
 
  implicit none

  private
  public :: compute_interface_flux_upwind
  public :: compute_flux_jacobian_upwind

  contains

  subroutine compute_interface_flux_upwind(q, dt, flux)
    
    real(r8), intent(in) :: q(1-nhalo:nx+nhalo)
    real(r8), intent(in) :: dt
    real(r8), intent(out) :: flux(0:nx)
    integer i

    select case (trim(space_scheme))
    case ('upwind1')
      do i = 0, nx
        flux(i) = merge(u * q(i), u * q(i+1), u > 0.0_r8)
      end do
    case ('upwind3')
      do i = 0, nx
        flux(i) = u / 12.0_r8 * (7.0_r8 * (q(i+1) + q(i)) - (q(i+2) + q(i-1))) - &
             abs(u) / 12.0_r8 * (3.0_r8 * (q(i+1) - q(i)) - (q(i+2) - q(i-1)))
      end do
    case default
      print*, 'Unknown space_scheme: ', trim(space_scheme)
      stop
    end select
    
  end subroutine compute_interface_flux_upwind

  subroutine compute_flux_jacobian_upwind(q, dt, jac)

    real(r8), intent(in) :: q(1-nhalo:nx+nhalo)
    real(r8), intent(in) :: dt
    real(r8), intent(out) :: jac(1:nx, 1:nx)
    integer :: i, idx_m1, idx_p1, idx_m2, idx_p2

    jac = 0.0_r8
    select case (trim(space_scheme))
    case ('upwind1')
      if (u > 0.0_r8) then
        do i = 1, nx
          idx_m1 = pbc_idx(i-1, nx)
          jac(i, idx_m1) = -u
          jac(i, i     ) = u
        end do
      else
        do i = 1, nx
          idx_p1 = pbc_idx(i+1, nx)
          jac(i, i     ) = -u
          jac(i, idx_p1) = u
        end do
      end if
    case ('upwind3') ! Wicker and Skamarock, 2002
      do i = 1, nx
        idx_m2 = pbc_idx(i-2, nx)
        idx_m1 = pbc_idx(i-1, nx)
        idx_p1 = pbc_idx(i+1, nx)
        idx_p2 = pbc_idx(i+2, nx)
        jac(i, idx_m2) = -(-u / 12.0_r8 - abs(u) / 12.0_r8)
        jac(i, idx_m1) = -u / 12.0_r8 - abs(u) / 12.0_r8 - &
                        (7.0_r8*u / 12.0_r8 + 3.0_r8 * abs(u) / 12.0_r8)
        jac(i, i     ) = 7.0_r8*u / 12.0_r8 + 3.0_r8 * abs(u) / 12.0_r8 - &
                        (7.0_r8*u / 12.0_r8 - 3.0_r8 * abs(u) / 12.0_r8)
        jac(i, idx_p1) = 7.0_r8 * u / 12.0_r8 - 3.0_r8 * abs(u) / 12.0_r8 -&
                         (-u / 12.0_r8 + abs(u) / 12.0_r8)
        jac(i, idx_p2) = -u / 12.0_r8 + abs(u) / 12.0_r8
      end do
    case default
      print*, 'Unknown space_scheme: ', trim(space_scheme)
      stop
    end select
    jac = -jac / dx

  end subroutine compute_flux_jacobian_upwind

end module upwind_mod